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Abstract 

In this paper we consider a heterogeneous network of Haken-Kelso- 
Bunz (HKB) nonlinear oscillators coupled through both linear and 
nonlinear interaction protocols. While some work exists on a system 
made up of only two nonlinearly coupled HKB oscillators as a model 
of human dynamics during interpersonal coordination tasks, the prob¬ 
lem of considering a network of three or more HKBs has not been 
fully investigated. The aim of our work is to study convergence and 
synchronization in networks of HKB oscillators as a paradigm of coor¬ 
dination in multiplayer games. Convergence results are obtained under 
the assumption that the network is connected, simple and undirected. 
Analytical results are obtained to prove convergence when the oscilla¬ 
tors are coupled diffusively. All theoretical results are illustrated via 
numerical examples. Finally, the effects of adding an external entrain¬ 
ment signal to all the agents in the network are analyzed and a model 
to account for them is proposed. 


1 Introduction 

Interpersonal coordination and synchronization between the motion of two 
individuals have been extensively studied over the past few decades [iniiia 

EDIS]- 

Synergetic movements of two or more people mirroring each other fre¬ 
quently occur in many activities such as handling objects, manipulating a 
common workpiece, dancing, choir singing and movement therapy miEiiEa 
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[Mll33| . It is of great importance to reveal not only the effects of mirroring 
movements among people on human physiological and mental functions, but 
also to deeply understand the link between intrapersonal and interpersonal 
coordination. In social psychology, it has been shown that people prefer to 
team up with others possessing similar morphological and behavioral fea¬ 
tures, and that they tend to coordinate their movement unconsciously mm- 
Moreover, much evidence suggests that motor processes caused by interper¬ 
sonal coordination are strictly related to mental connectedness. In partic¬ 
ular, motor coordination between two human subjects contributes to so¬ 
cial attachment particularly when the kinematic features of their movement 
share similar patterns [32ll36j . 

In order to explain the experimental observations of human interper¬ 
sonal coordination, mathematical models are usually derived to capture the 
key features of the observed behavior. A classical example is the so-called 
HKB oscillator which was introduced in |10] to explain the transition from 
phase to antiphase synchronization in bimanual coordination experiments 
(for more details see |13lll7] i. The HKB nonlinear oscillator was shown to 
be able to capture many features of human coordination even beyond the 
bimanual synchronization experiments it was derived to explain. For ex¬ 
ample, HKB oscillators were used in [391140] in the context of the mirror 
game [23], often presented as an important paradigmatic case of study to 
investigate the onset of social motor coordination between two players imi¬ 
tating each other’s hand movements. Furthermore, in m the authors take 
inspiration from the dynamic clamp of cellular and computational neuro¬ 
science in order to probe essential properties of human social coordination 
by reciprocally coupling human subjects to a computationally implemented 
model of themselves (HKB oscillator), referred to as Virtual Player (VP). 
Such concept, namely the human dynamic clamp (HDC), was further in¬ 
vestigated and developed in |1] in order to cover a broader repertoire of 
human behavior, including rhythmic and discrete movements, adaptation to 
changes of pacing, and behavioral skill learning as specified by the VP. Be¬ 
sides, HKB oscillators were also used in m in order to capture the rhythmic 
coordination between two individuals swinging hand-held pendulums, in [34] 
in order to model spontaneous interpersonal postural coordination between 
two human people and account for the competition between the coupling to 
a visual target to track and the coupling to the partner, in |26| in order to 
qualitatively explain interpersonal movement synchronization between two 
human beings involved in rhythmic paradigms, and in [T] in order to account 
for the frequency detuning of the phase entrainment dynamics of two people 
involved in interlimb coordination. 

While coordination of two human players has been studied in numer¬ 
ous previous investigations, the case of multiple human players has been 
seldom studied in the existing literature, due to a combination of practical 
problems in running the experiments and lack of a formal method able not 


2 



only to model the considered scenario but also to quantify and character¬ 
ize the synchronization level of the ensemble. Multiplayer games involve a 
group of three or more people engaged in a communal coordination task. 
The variety of scenarios that can be considered is vast due to the countless 
activities the players might be involved in (limb movements, finger move¬ 
ments, head movements, walking in a crowd, or more in general music and 
sport activities), the many ways in which the participants can interact and 
communicate with each other and the different ways all the players can be 
physically located with respect to each other while performing the specified 
task. 

Some of the existing works on coordination of multiple human players in¬ 
clude studies on choir singers during a concert m , rhythmic activities as for 
example ’’the cup game” and marching tasks [12], rocking chairs |6l|25] and 
coordination of rowers’ movements during a race m- In these papers the 
authors provide several experimental results in order to analyze the behavior 
of a group of people performing some coordinated activities, but a rigorous 
mathematical model capable of capturing the observed results and explain¬ 
ing the features of the movement coordination among them is still missing. 
In particular, in |6| the authors study both unintentional and intentional 
coordination by asking the players to try and synchronize the oscillations of 
the rocking chairs with their eyes shut or open. Synchronization is observed 
to spontaneously emerge when players observe each other’s movements. An¬ 
other study in which multiplayer activities are analyzed but a mathematical 
model is missing is carried out in [38|: the authors use the symmetric Hopf 
bifurcation theory, which is a model-independent approach based on coupled 
oscillators, to investigate the synchronized patterns of three people during 
sport activities. 

Further results about multiplayer activities deal with spontaneous group 
synchronization of arm movements and respiratory rhythms. For example, 
in [2] the authors test whether pre-assigned arm movements performed in 
a group setting spontaneously synchronize and whether synchronization ex¬ 
tends to heart and respiratory rhythms. In their study no explicit directions 
are given on whether or how the arm swingings are to be synchronized among 
participants, and experiments are repeated with and without external cues. 
Interestingly, when an external auditory rhythm is present, both motor and 
respiratory synchronization is found to be enhanced among the group. Also 
the overall coordination level is observed to increase when compared to that 
detected when the same experiments are again carried out in its absence. 

While in [2] no mathematical model is presented, in [9] the effects of 
an external visual signal, when considering the postural sway of a human 
being, is explicitly modeled by introducing a linear oscillator with adaptive 
parameters depending on the frequency of the exogenous visual signal itself. 
Moreover, when considering the social postural coordination between two 
human beings |34j . another model for an external visual stimulus is presented 
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by adding a unidirectional coupling term to the dynamics of a nonlinear 
oscillator. 

The main objective of this paper is to propose and analyze a model able 
to account for the onset of movement synchronization in multiplayer scenar¬ 
ios and explain some of the features observed experimentally in the existing 
literature. Specifically, we consider a heterogeneous network of HKB nonlin¬ 
ear oscillators as a good model of multiplayer coordination and, as already 
done in [22] for the case of two agents only, we regard it as a synchronization 
problem. Each equation is used to model the movement of a different player 
and is therefore characterized by a different set of parameters to account for 
human-to-human variability. The effect of different interaction models, lin¬ 
ear and nonlinear, will be investigated to understand under what conditions 
synchronization is observed to emerge. Our analysis suggests that bounded 
synchronization is indeed a common emergent property in these networks 
whose occurrence can also be accounted for analytically in a number of dif¬ 
ferent scenarios. Also, as expected from existing theoretical results, we find 
that the structure of the interactions among players has an effect on the 
coordination level detected in the network. 

Furthermore, the effects of adding an external sinusoidal signal will be 
studied in order to understand whether synchronization can be improved by 
means of an entrainment signal [28|. Our analysis suggests that the syn¬ 
chronization level of the ensemble can indeed increase when the oscillation 
frequency of the external signal is similar to the natural angular velocity of 
the agents in the network. However, in all the other cases, the external sig¬ 
nal acts as a disturbance and leads to a decrease of the coordination among 
the agents. 

We wish to emphasize that the study reported in this paper will form 
the basis of future experimental investigations which are currently being 
planned. 

The paper is organized as follows. In Sect. [2] some notation that shall be 
used in later sections is introduced. In Sect. Ejthe equation that describes 
the network is presented, in terms of both internal dynamics of each agent 
and coupling protocol thanks to which they can interact with each other. 
In Sect. Hlsome metrics are introduced to characterize the quality and the 
level of coordination in human groups. In Sect. [5] a testbed scenario of 
multiplayer coordination in networks of human people is presented, while 
in Sect. 15.11 the key synchronization features experimentally observed are 
reproduced by considering a heterogeneous network of HKB oscillators, and 
the effects of three different coupling strategies thanks to which they are 
interconnected are explored. In Sect. O the effects of adding an external 
entrainment signal is analyzed with respect to the overall synchronization 
level of the network. In Sect. [7] global bounded synchronization of the 
network when its nodes are connected through a linear diffusive coupling 
protocol is analytically proven to be achieved, and some numerical examples 
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are provided in order to both illustrate the effectiveness of our analysis and 
to show that bounded synchronization can be achieved also when considering 
different couplings. Finally, in Sect. [8] a summary of our results and some 
possible future developments are presented. 

2 Preliminaries and background 

We denote with 0 the Kronecker product between two matrices. The oper¬ 
ator Afc (•) defined over a matrix indicates the k-th. eigenvalue of the matrix 
itself, and Am (•) indicates its maximum eigenvalue when the matrix is real 
and symmetric and as a consequence all the eigenvalue are real as well. 

A graph is a tuple Q = {V, Al} defined by a set of nodes V = {1,..., N} 
and a set of edges Al C V x V. A graph is said to be undirected if (i,j) G 
A (j, i) G Al. In an undirected graph, two nodes i and j are said to 

be neighbors if (f,j) G Al. The matrix A = {aij} G where 

> 0, if {i,j) a-i'e neighbors 
= 0, otherwise 

is called adjacency matrix, and aij > 0 is called strength of the interaction 
between the pair (i,j). In particular, a graph is said to be unweighted 
if the interaction between two neighbors is equal to 1. A path between 
nodes h and A; is a sequence of nodes, with h and k as endpoints, such that 
every two consecutive nodes are neighbors. A graph is said to be simple if 
an = 0 Vf G V, while it is said to be connected if there exists a path between 
any two of its nodes. The matrix L = G defined as 

Ylk=i O'ik, 'Ai = j 

-Oij, 'AiAj 

is called Laplacian matrix of the graph (or simply Laplacian). The Lapla- 
cian of any simple undirected graph is symmetric with zero row sum and 
is a positive semidefinite matrix with as many null eigenvalues as there are 
components in the graph. In particular, a connected graph has only one null 
eigenvalue. 

Throughout the paper we shall consider a connected simple undirected 
network of N agents assuming that any two players interact symmetrically 
with one another. 

Before analyzing a multiplayer scenario, it is worth considering the sim¬ 
pler case of only two human players interacting with each other. The system 
that can be used to model the interaction between them is described as fol¬ 
lows [ZIES]: 
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(2) 


! xi + [ax\ + (5x\ - 7 ) xi +ujIxi = I{xi,X 2 ) 

X2 + {axl + ( 3 x 2 -1) X2+ ^2^2 = I{x2,Xi) 

where Xj G M denotes the position of the i-th player, with i = 1,2. The 
right-hand side of both equations represents the coupling term between the 
two players: in particular 

I {w, z) := [a + b{w — z)‘^] {w — z) (3) 

The term (ax? + /3x? — 7 ) Xi represents the nonlinear damping of the os¬ 
cillatory movement of player i. Specifically, the sign of 7 determines whether, 
in the absence of coupling, the oscillation is persistent (7 > 0 ) or vanishes 
(vice versa) as time goes by: it is trivial to verify this by studying the stabil¬ 
ity of the origin and checking the sign of the eigenvalues of the Jacobian of 
the system. Moreover, a and (3 determine the amplitude of such oscillation, 
while LJi is related to its frequency. It has been proven that this model of 
two nonlinearly coupled oscillators account for the features observed during 
experimental data in bimanual experiments (see [lOj for further details). 


3 Human to human coordination as a synchroniza¬ 
tion problem 

In the introduction of this paper we have pointed out that the dynamics 
of two coupled HKB oscillators has been used to describe different kinds 
of interpersonal coordination tasks between two people, including bimanual 
coordination experiments, mirror game, social postural coordination and 
rocking chairs. According to the particular scenario considered, the state 
vector of each oscillator is used to represent position and velocity of the 
particular body part of interest of either of the players (finger, hand, head, 
and so forth). Following the same approach, we can consider a scenario in 
which more than two human beings are performing a multiplayer coordina¬ 
tion task, as for example arm or hand rhythmic movements, rocking chairs, 
head tracking of a visual target and so on. In these cases, the state vector 
of each node represents position and velocity of the particular body part of 
interest of each player. Therefore, the dynamics of each player when moving 
in isolation will be described by an HKB equation: 


fi{t,Xi) 


Xi2 

-{aix\ -h Axf^ - li)xi2 - w^xq 


(4) 


where x, = [xq Xjj]'^ G is the state vector, with xq,Xj 2 representing 
position and velocity of the i-th human player, respectively. 

To model the interaction between different players we assume that the 
dynamics of each of them is affected by some coupling function Ui which 


6 




depends on the state of its neighbors. In what follows we will explore the 
effects of three possible selections for such a function. In particular, we are 
interested in analyzing the differences of the results provided by all of them 
and understanding which one leads to synchronization features which are 
the closest to the ones observed in some existing work about group syn¬ 
chronization of networks of several human people involved in a coordination 


task |25j . 


1. Full state coupling. With this kind of coupling, we assume that play¬ 
ers adjust both their velocities and accelerations proportionally to the 
average mismatch between theirs and those of their neighbors. Math¬ 
ematically, we have: 


N 



In particular, A/i > 0 is the number of neighbors of node i, while c > 0 
is the coupling strength among the agents. 

2. Partial state coupling. Next, we explore the case where players only 
adjust their accelerations according to the position and velocity mis¬ 
matches from their neighbors: 



In particular, A/i > 0 is the number of neighbors of node i, while 
ci,C 2 > 0 represent the position and the velocity coupling strengths, 
respectively. 

3. HKB coupling. Finally we consider an interaction model which is the 
direct extension to multiplayer coordination problems of the interac¬ 
tion function used in the classical HKB set up to model coordination 
between two players mm- Specifically we choose the following non¬ 
linear function: 


M aij[a + b{xi^ - - xj^) 


Once again, A/i > 0 is the number of neighbors of node i, while c > 0 
represents the coupling strength among the agents. 

The resulting network model describing the interaction of a group of N 
players can then be written as 
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( 8 ) 


Xi{t) 


-{aiX 


+ 


^12 

- 7i)xi2 - uj'jxi 


+ Ui{t) G 


where the coupling function Ui can be chosen as one of those listed above. We 
now explore under what conditions coordination, and hence synchronization, 
emerges for each of the three scenarios of interest. 

We wish to emphasize that, since the node parameters are heterogeneous, 
complete synchronization as defined in [20] cannot be achieved. We will 
consider instead the case where bounded synchronization, as defined below, 
emerges. Namely, we define the average trajectory as 


and the tracking error as 


x{t) : = 


i=i 


(9) 


ei{t) := Xi{t) — x{t) yt > 0 ,i = 1 , N (10) 

We also define the parameters vector for each node i as := [a^ Pi ji uji]'^ G 
and we introduce the stack vectors x{t) := [xi{t)'^ X 2 {t)'^ ... x^it)'^]'^ G 
and e{t) := ■■■ eAr(i)'^]^ £ and the error norm 

rj(t) := ||e(f)|| G M,Vf > 0, where || • || indicates the Euclidean norm. 

We say that a network of HKB oscillators achieves coordination if and 
only if 


lim r]{t) < e (11) 

t^OO 

for any initial condition xt^ and parameter vector -dj of the nodes in the 
network, where e > 0 is a sufficiently small constant. 

4 Synchronization metrics 

In order to quantify and analyze the synchronization level in a network of 
more than two agents, we use the metrics introduced in [25] to characterize 
the quality and the level of coordination in human groups. 

Let Xk{t) G M Vf G [0,T] be the continuous time series representing the 
motion of each agent, with k G [1, iV], where N is the number of individuals 
and T is the duration of the experiment. Let Xk{ti) G M, with k G [1, N] and 
i G [1, Nt], be the respective discrete time series of the k-th agent, obtained 
after sampling Xk{t), where Nt is the number of time steps and AT := ^ 
is the sampling period. Let 0k{t) G [—vr,7r] be the phase of the k-ih. agent, 
which can be estimate by making use of the Hilbert transform of the signal 
Xk{t) [l8|. We define the cluster phase or Kuramoto order parameter, both 
in its complex form q'{t) G C and in its real form q{t) G [—7r,7r] as 
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^ k=l 

q{t) := atan2 {^{q'{t)),^{q'{t))) (13) 

which can be regarded as the average phase of the group at time t. 

Let 4>k{t) ■= 0k{t) — q{t) be the relative phase between the /c-th par¬ 
ticipant and the group phase at time t. We can define the relative phase 
between the /c-th participant and the group averaged over the time interval 
both in its complex form G C and in its real form 
as 


^'k--= 


1 

T 




Nt 

_L p{j4>k{U)} 

Nt ^ 

1=1 


(14) 


:= atan2 (15) 

In order to quantify the degree of synchronization for the /c-th agent 
within the group we define the following parameter 


Pk:=\^'k\ G[0,1] (16) 

which simply gives information on how much the fc-th agent is synchronized 
with the average trend of the group. The closer pk is to 1, the better the 
synchronization of the /c-th agent itself. 

In order to quantify the synchronization level of the entire group at time 
t we define the following parter 


Pg{t) ■ = 


1 

iv 


N 

V p{Mk(i)-'t>k]} 


k=l 


G [0,1] 


(17) 


which simply represents the group synchronization: the closer Pg{t) is to 1, 
the better the synchronization level of the group at time t. Its value can be 
averaged over the whole time interval [0, T] in order to have an estimate of 
the mean synchronization level of the group during the total duration of the 
performance; 


pg:= ^ pg{t) ^ ^ PgiU) G [0,1] (18) 

Finally if we denote with (kd^yif) '■= 0k{t) — 0k'{t) the relative phase 
between two participants in the group at time t, it is possible to estimate 
their dyadic synchronization, that is the synchronization level between par¬ 
ticipants k and k' over the whole round: 
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G[0,1] (19) 


Nt 

Nt ^ 

1=1 

It is worth pointing out that high dyadic synchronization levels can co¬ 
exist with low group synchronization values. 

5 Testbed example 

As a testbed scenario we consider the synchronization of rocking chairs mo¬ 
tion studied in [25]. In particular, participants sit on six identical wooden 
rocking chairs disposed as a circle and are supposed to rock them in two 
different conditions: 

1. Eyes closed: participants are required to rock at their own preferred 
frequency while keeping their eyes closed; 

2. Eyes open: participants are required to rock at their own preferred 
frequency while trying to synchronize their rocking chair movements 
as a group. 

In the eyes closed condition the participants are not visually coupled, 
meaning that the oscillation frequency of each of them is not influenced by 
the motion of the others, whilst in the eyes open condition each player is 
asked to look at the middle of the circle in order to try and synchronize their 
motion with the one of the others. The six participants first perform a trial 
while keeping their eyes closed, then perform two eyes open trials, namely 
T1 and T2; each of the three trials lasts 3 minutes. 

In Fig. [T] some results about the typical trend of the group synchroniza¬ 
tion Pg{t) and its mean value and standard deviation are represented for 
each of the three aforementioned trials. In particular, in Fig. [T^ we can 
observe that the mean value pg of the group synchronization, represented as 
a circle, is around 0.4 in the eyes closed condition, while it is around 0.85 
in the eyes open condition: this means that, when the participants are not 
visually coupled, as expected synchronization does not emerge, whilst when 
visually coupled and explicitly told to rock their own chair movements as a 
group, the synchronization level significantly increases. In Fig. [Tja we can 
observe that in the eyes closed condition the amplitude of the oscillations of 
the group synchronization is higher than the one obtained in the eyes open 
condition. 

In Table [U we show typical values of the degree of synchronization pk 
of the participants involved in the rocking chairs experiments, both for the 
eyes closed and the eyes open condition. It is easy to see that, as expected, 
the value of p^ is much higher for almost all the participants when they 


Pd.k,k' ■ 




dt 
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(a) mean value and standard deviation (b) typical trend 

Figure 1 Group synchronization in the rocking chairs experiments of [25] - 
T1 and T2 refer to two different trials of the eyes open condition 


are visually coupled. Interestingly enough, agent 6 does not undergo an 
improvement of pQ with respect to the eyes closed condition, meaning that 
such participant has more trouble synchronizing with the group compared 
to the other ones. 

Table 1 Degree of synchronization of the participants in the rocking chairs 
experiments of [25| - EC: eyes closed, EO: eyes open 


Participant 

EC {pg = 0.36) 

EO {pg = 0.80) 

1 

0.36 

0.95 

2 

0.34 

0.92 

3 

0.30 

0.95 

4 

0.35 

0.88 

5 

0.34 

0.67 

6 

0.40 

0.37 


In Table [2] we show typical values of the dyadic synchronization of 
the participants involved in the rocking chairs experiments for the eyes open 
condition. As expected, the lowest values are the ones obtained with respect 
to the participant that most struggled to synchronize with the group, that 
is participant 6. 

5.1 Modeling results 

In this section we uncover the synchronization features that the three differ¬ 
ent interaction protocols introduced earlier lead to, with respect to the rock¬ 
ing chairs experiments introduced earlier as a testbed scenario |25] . We will 
explore whether and how the model of coupled HKB oscillators we propose 
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Table 2 Dyadic synchronization of the participants in the rocking chairs 
experiments of [25] for the eyes open condition [pg = 0.80) 


Participants 

2 

3 

4 

5 

6 

1 

0.87 

0.86 

0.81 

0.63 

0.19 

2 

— 

0.85 

0.78 

0.59 

0.21 

3 

— 

— 

0.82 

0.61 

0.21 

4 

— 

— 

— 

0.50 

0.18 

5 

— 

— 

- 

- 

0.14 


in this paper can reproduce the key features of the observed experimental 
results. In so doing we will explore: 

• the effects of choosing different coupling functions; 

• how varying the coupling strength affects the synchronization level of 
the agents. 

In what follows we simulate a heterogeneous network of = 6 HKB 
oscillators whose parameters and initial values are heuristically set as de¬ 
scribed in Table [3] and we set T = 200s. We suppose that the network is 
simple, connected, unweighted and undirected and we assume that each node 
is connected to all the others (complete graph), which we believe well rep¬ 
resents the topology implemented in the rocking chairs experiments of [25] 
for the eyes open condition. 

Table 3 Numerical simulations - parameters and initial values for a network 
oi N = 6 HKB oscillators 


Nodes ai 

A 

li 

Ui 

Xi{0) 


1 

0.46 

1.16 

0.58 

0.31 

[-1.4, 

+0.3] 

2 

0.37 

1.20 

1.84 

0.52 

[+1.0, 

+0.2] 

3 

0.34 

1.73 

0.62 

0.37 

[-1.8, 

-0.3] 

4 

0.17 

0.31 

1.86 

0.41 

[+0.2, 

-0.2] 

5 

0.76 

0.76 

1.40 

0.85 

[+1.5, 

+0.1] 

6 

0.25 

0.86 

0.56 

0.62 

[-0.8, 

-0.1] 


In particular, since we are interested in replicating the key features of 
the rocking chairs experiments for both the conditions (eyes open and eyes 
closed), in Fig. [2] we show the group synchronization obtained with and 
without interaction protocol. In particular, in Fig. [2^ we show mean value 
and standard deviation of Pg{t): in each column, they are shown both for all 
the three interaction protocols presented earlier and for the case in which 
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(a) mean value and standard deviation (b) typical trend 

Figure 2 Group synchronization in a heterogeneous unweighted complete 
graph of = 6 HKB oscillators - a NC: no coupling, FSC: full state 
coupling (c = 0.15), PSC: partial state coupling (ci = C 2 = 0.15), HKB: 
HKB coupling {a = b = —l,c = 0.15) - b black dashed line: no coupling, 
red solid line: full state coupling 


the nodes are not connected, respectively. The mean value is indicated with 
a circle, while the standard deviation from it is indicated with a vertical 
bar whose extremities are delimited by two horizontal lines. In particular, 
if we denote with pg the mean value of the the group synchronization Pg{t) 
obtained as defined in Eq. [18] for each of the four aforementioned cases after 
a simulation of duration T and with its standard deviation, respectively, 
we have that: 


= 




E 


(pgiti)- Pgf 


( 20 ) 


It is easy to observe that, in absence of connections among the nodes, 
which corresponds to = 0 Vi G [1;A^]; the group synchronization has 
a mean value approximately equal to 0.4, while it significantly increases 
(approximately 0.9) when connecting the nodes with any of the three inter¬ 
action protocols introduced above. These results confirm the observations 
previously made for a network of six human people involved in rocking chairs 
experiments (see Fig. [1^). In particular, we have chosen c = 0.15 for the full 
state interaction protocol, ci = C 2 = 0.15 for the partial state interaction 
protocol, and a = b = —1, c = 0.15 for the HKB interaction protocol. 

In Fig. [2}3 we show the time evolution of the group synchronization Pg{t) 
when the nodes are not connected at all (black dashed line) and when they 
are connected through a full state interaction protocol (red solid line): for 
the sake of clarity we do not show the trend of pg (t) obtained with a partial 
state and an HKB interaction protocol since they are qualitatively analogous 
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to the one obtained with a full state interaction protocol. Our simulation 
results are able to reproduce another key feature observed in [25] : when the 
nodes are uncoupled, which corresponds to the eyes closed condition, the 
amplitude of the oscillations of Pg{t) is higher than the one obtained when 
the nodes are coupled instead, which corresponds to the eyes open condition 
(see Fig. Eb). 

Then, in Table 0] we show the degree of synchronization obtained for 
each node of the network, both in absence of coupling among the agents 
and in its presence. It is easy to see that for each node k in the network, 
Pk has much higher values when any of the three interaction protocols is 
introduced, confirming what observed in [25] when asking the participants to 
rock their chairs while keeping their eyes open rather than closed. Moreover, 
we are able to reproduce another interesting feature: despite the group 
synchronization assuming high values when the human players are visually 
coupled, there might be some agents that struggle to keep up with the 
general trend of the group, therefore showing lower values in terms of pk 
(node 5 in our simulations). 

Table 4 Degree of synchronization of the nodes in a heterogeneous un¬ 
weighted complete graph of = 6 HKB oscillators - no coupling (NC), 
full state coupling (FSC) with c = 0.15, partial state coupling (PSC) with 
Cl = C 2 = 0.15 and HKB coupling (HKB) with a = b = —1, c = 0.15 


Node 

NC 

FSC 

PSC 

HKB 

1 

0.42 

0.95 

0.93 

0.97 

2 

0.38 

0.92 

0.94 

0.96 

3 

0.45 

0.98 

0.95 

0.97 

4 

0.41 

0.98 

0.96 

0.98 

5 

0.33 

0.33 

0.33 

0.49 

6 

0.36 

0.98 

0.97 

0.98 


Furthermore, in Table 0 we show the dyadic synchronization p^^ for 
all the possible couples of nodes in the network: again, our simulation re¬ 
sults confirm what observed for the rocking chairs experiments. Indeed, the 
couples of nodes with lower dyadic synchronization correspond to pairs in 
which at least either of the two nodes had trouble synchronizing with the 
general trend of the group (node 5 in our simulations). For the sake of 
clarity we show p^^ only when connecting the nodes through a full state 
interaction protocol, since analogous results can be obtained also for the 
other two strategies introduced earlier in this paper. 

It is easy to foresee that, regardless of the interaction protocol the nodes 
are connected through, the value of the coupling strength has a direct impact 
on the group synchronization in terms of its mean value and its standard 
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Table 5 Dyadic synchronization in a heterogeneous unweighted complete 
graph of = 6 HKB oscillators - full state coupling (c = 0.15) 


Nodes 

2 

3 

4 

5 

6 

1 

0.91 

0.98 

0.94 

0.38 

0.98 

2 

— 

0.92 

0.96 

0.43 

0.93 

3 

— 

— 

0.97 

0.39 

0.99 

4 

— 

— 

— 

0.40 

0.98 

5 

— 

- 

- 

— 

0.41 




(a) group synchronization (b) group synchronization - zoom 

Figure 3 Mean value and standard deviation of the group synchronization 
in a heterogeneous unweighted complete graph oi N = 6 HKB oscillators 
for different values of the coupling strength c - full state coupling 


deviation. We now show how quantitatively pg varies as the the coupling 
strength varies for all the three interaction protocols introduced earlier in 
this paper, when considering a heterogeneous unweighted complete graph of 
N = 6 HKB oscillators whose parameters and initial values are described in 
Table [3l Moreover, we set T = 200s. 

From Fig. [3] to 0 we show mean value and standard deviation of the 
group synchronization Pg{t) obtained for different values of the coupling 
strength when considering full state coupling, partial state coupling and 
HKB coupling as interaction protocols, respectively. In particular, the blue 
solid line refers to the mean value of pg (t) , while the red dashed lines indicate 
the variation or dispersion from it. 

From Fig. [3^ it is clear that, when considering a full state coupling 
as interaction protocol, the group synchronization increases as the coupling 
strength c increases: in particular, in order for the network to well syn¬ 
chronize it is sufficient a relatively small value for the coupling strength 
(c ~ 0.15, see Fig. [ 8 ) 3 ). In terms of multiplayer games, this means that 
the stronger the influence that each player has on the others, the better the 
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(a) C2 = 0, Cl variable (b) Cl = 0, C2 variable 

Figure 4 Mean value and standard deviation of the group synchronization 
in a heterogeneous unweighted complete graph of = 6 HKB oscillators for 
different values of the coupling strengths ci and C 2 - partial state coupling 


overall synchronization of the human participants. 

In Fig. [3^ we show how, when considering a partial state coupling as 
interaction protocol, the group synchronization varies for increasing values 
of Cl while keeping C 2 constantly equal to 0, and vice versa in Fig. 0 ) 3 . As we 
can see, the influence that C 2 has on the group synchronization is stronger 
than the one provided by ci. In terms of multiplayer games, this means 
that human players react better to changes in the velocity of their neighbors 
rather than in their position. This results is confirmed also in Fig. [ 6 ^ in 
which we show how the mean value of the group synchronization changes as 
Cl and C 2 are simultaneously varied (darker colors refer to lower values of the 
average group synchronization, whilst lighter ones refer to higher values). 

Finally from Fig. it is clear that, when considering an HKB coupling 
as interaction protocol while keeping a and h constantly equal to — 1 , the 
group synchronization increases as the coupling strength c increases. In 
particular, like in the case of a full state interaction protocol, in order for the 
network to well synchronize it is sufficient to choose a relatively small value 
for the coupling strength (c ~ 0.15, see Fig. [Sh). In terms of multiplayer 
games, this means that the stronger the influence that each player has on 
the others, the better the overall synchronization of the human participants. 
This results is confirmed also in Fig. [ 6)3 in which we show how the mean 
value of the group synchronization changes as a and h are simultaneously 
varied while keeping c constantly equal to 1 (darker colors refer to lower 
values of the average group synchronization, whilst lighter ones refer to 
higher values). As we can see, as the values of |a| and |5| increase, then the 
average of the group synchronization increases as well. 
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(a) group synchronization (b) group synchronization - zoom 


Figure 5 Mean value and standard deviation of the group synchronization 
in a heterogeneous unweighted complete graph oi N = 6 HKB oscillators 
for different values of the coupling strength c while keeping a = b = —1 
constant - HKB coupling 



0 0.2 0.4 0.6 0.8 0 0.1 0.2 03 0.4 03 0.6 

Cl -a 


(a) partial state coupling - ci,C2 variable (b) HKB coupling - a, b variable while keep¬ 
ing c = 1 constant 

Figure 6 Mean value of the group synchronization Pg{t) in a heterogeneous 
unweighted complete graph of = 6 HKB oscillators for different values of 
the coupling strengths 
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6 Entrainment of the network 


In this section we analyze the effects on the group synchronization of the 
network of adding an external sinusoidal signal to the dynamics of each 
node. Our main objective is understanding whether, and possibly under 
what conditions, such entrainment signal leads to a better synchronization 
level of a heterogeneous network of HKB oscillators with respect to the case 
in which the signal is absent. This will help us understand whether an 
external auditory or visual stimulus can improve the coordination level in 
multiplayer games when considering networks of human people involved in 
some synchronization task. 

Following the approach of [28] , we model such a scenario in the following 
way: 




-{aix: 


+ A4 


- -ii)Xi^ - UJ-X 




+ c 


+ Ui 


( 21 ) 


where ({t) = ^^sin(a;(^f) represents the entrainment signal and Ui{t) one of 
the interaction protocols introduced earlier in this paper. 

We introduce the entrainment index pE € [0,1] in order to quantify the 
overall synchronization level between the network and the external signal 

City. 


PEk ■ = 


1 r e^[ek{t)-8c{t)] 

T Jo 


PE '■ = 


1 

iv 


N 

^PEk 


k=l 


( 22 ) 


where 6kit) is the phase of the /c-th node, Of^it) is the phase of Cit), T is the 
duration of the experiment and N is the number of nodes in the network. 
The closer pE is to 1, the better the synchronization of the group with the 
entrainment signal. 

In what follows we simulate a heterogeneous network of = 6 HKB os¬ 
cillators whose parameters and initial values are heuristically set as described 
in Table [3] and we set T = 200s. We suppose that the network is simple, 
connected, unweighted and undirected and we assume that each node is con¬ 
nected to all the others (complete graph), which we believe well represents 
the topology implemented in the rocking chairs experiments of [25]. 

In Fig. [7| we show the entrainment index for different values of the 
frequency and the amplitude of the entrainment signal Cit) when 
considering a full state coupling as interaction protocol with c = 0.15 (darker 
colors refer to lower values of pE, whilst lighter ones refer to higher values). 
It is easy to see that, for each value of the entrainment index increases 
as A(^ increases as well, meaning that the network better synchronizes with 
Cit) for increasing values of its amplitude. Moreover, for a given value of 
the highest values of pE are achieved when the frequency of the entrainment 
signal is close to the average value H of the natural frequencies uji of the 
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Figure 7 Entrainment index in a heterogeneous unweighted complete graph 
of A'" = 6 HKB oscillators - full state coupling (c = 0.15) 


nodes (in this case Q ~ 0.5 ). These results conhrm the findings of |30[I35| . 
in which it is shown that spontaneous unintentional synchronization between 
the oscillation of a handheld pendulum swung by an individual and of an 
external sinusoidal stimulus (which corresponds to our external entrainment 
signal) emerges only when the frequency of the signal itself is similar to the 
preferred frequency of the player. For the sake of brevity we do not show 
how pe varies as ujq and Aq vary as well when considering partial state 
coupling and HKB coupling as interaction protocols, since we obtain results 
which are analogous to the ones shown in Fig. [7] for a full state coupling. 

In Fig. [8] we show mean value and standard deviation of the group syn¬ 
chronization Pg{t) when considering different parameters of the entrainment 
signal (green line: no entrainment signal, red line: u)(^ = 0.1,74^ = 0.1, 
blue line: = 0.3, Aq = 0.2, black line: ujq = 0.5, ^4^^ = 0.3) for all 

the three coupling protocols we have presented {FSC: full state coupling 
(c = 0.15) - PSC: partial state coupling (ci = C 2 = 0.15) - HKB\ HKB 
coupling [a = b = —l,c = 0.15)). Since we are interested in understanding 
whether an additive external sinusoidal signal can improve the synchroniza¬ 
tion level of the network with respect to the case in which it is absent, the 
values of the coupling strengths chosen in these simulations for all the three 
interaction protocols are the same as the ones previously used in absence of 
entrainment signal (see Fig. [2^). From Fig. [8] it is easy to observe that, for 
all the three interaction protocols, the group synchronization of the network 
improves only when the entrainment index pE has high values (see black line 
compared to the green one). In the other two cases (blue line and red line), 
the entrainment signal acts as a disturbance for the dynamics of the nodes 
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1.1 



0 . 6 - 

O.sl-^^^- 

FSC PSC HKB 

Figure 8 Mean value and standard deviation of the group synchronization 
in a heterogeneous unweighted complete graph oi N = 6 HKB oscillators 
- FSC: full state coupling (c = 0.15) - PSC: partial state coupling (ci = 
C 2 = 0.15) - HKB: HKB coupling (a = b = —l,c = 0.15) - green line: no 
entrainment signal, red line: = 0.1, = 0.1, blue line: = 0.3, = 

0.2, black line: = 0.5, = 0.3 

and the group synchronization decreases. In terms of multiplayer games for 
networks of human people, this means that it is possible to further enhance 
the coordination level of participants only when the entrainment signal has 
an oscillation frequency which is close to the average of the natural oscilla¬ 
tion frequencies of the individuals involved and its amplitude is sufficiently 
high. 

7 Convergence analysis 

As anticipated earlier, other than finding a mathematical model able to re¬ 
produce features observed experimentally in some multiplayer games stud¬ 
ied in the existing literature, we are also interested in understanding under 
what conditions synchronization is observed to emerge. In particular, in 
this section we are going to show that global bounded synchronization can 
be analytically guaranteed for a heterogeneous network of diffusively cou¬ 
pled N HKB oscillators by making use of two different approaches, namely 
contraction theory and Lyapunov theory. 

7.1 Contraction theory 

Let I • I be a norm defined on a vector w G with induced matrix norm 11 • 11. 
As stated in [29], given a matrix P G the induced matrix measure is 
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defined as fJ,{P) '■= lim;j_^o+ 


(||/+feP||-i) 

h 


Definition 1. Let us consider the system w = F{t,w) definedMt > 0, Vw E 
C C M"". We say that sueh system is contracting with respect to a norm \ ■ \ 
with associated matrix measure /i(-) iff 


3 k > 0 : pL{J{w,t)) < —k, 'iw E C, Vt > 0 (23) 


where J is the Jacobian of the system. 

They key stage in the application of contraction theory to synchroniza¬ 
tion of networks of oscillators is the construction of the so-called virtual 
system [T3] . 

Definition 2. Let us consider a heterogenous network described by Eq. O 
The virtual system is defined as the system that has the trajectories of the 
nodes as particular solutions. 

Formally, the virtual system depends on the state variables of the os¬ 
cillators in the network and on some virtual state variables. Substitution 
of the state variables of a certain node i into the virtual ones returns the 
dynamics of the i-th node of the network itself. 

It is worth pointing out that virtual systems are originally defined for 
networks of identical systems in order to prove complete synchronization: 
indeed, it is possible to prove that if a virtual system is contracting, then 
limt_).oo r]{t) = 0. However, it is possible to define virtual systems also for 
networks of nonidentical oscillators by averaging the values of all the different 
parameters in order to prove bounded synchronization (see heterogeneous 
network of repressilators in m )• 

In [27] a simple algorithm is provided that allows to check whether the 
virtual system of a certain heterogeneous network of N agents is contract¬ 
ing, which leads to global bounded synchronization of the network itself. 
In particular, rather than verifying Eq. [231 in order to check whether the 
virtual system is contracting, the algorithm consists in checking the truth of 
some statements regarding the single elements of the Jacobian of the virtual 
system and imposing some conditions; 


1. build the Jacobian J of the virtual system; 

2. check whether the following statements are true or false 

• SI: J(i, i) < 0; 

• S2: J{i,i) = —pi, 0 < Pi < oo; 

• S3: J{i,j) / 0 ^ J(j,f) = 0; 

3. generate a set of conditions for synchronization (CFS) according to 
the truth or the falsity of the previous statements. 
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In particular, denoting with uq. the number of 0 elements in the i — th 
row of the Jacobian of the virtual system, the CFS are generated in the 
following way: 

. 51,52,53^ I J(i,j) I 

. 51,52,53^ I J(i,j)I 

. 51,52,53^ I I 

. 51,52,53^ I J(z,j)I 

Note that if statement 
system to be contracting. 

Theorem 1. Suppose to have a heterogeneous network of N HKB oscillators 
intereonneeted via a full state coupling as described in Eq. O Let us also 
assume that the network topology is a connected, simple, undirected and 
unweighted complete graph. If the following hypothesis is satisfied 

N-1, 9 ^ N-1 

—^ + w + 7 ) < c < —^ (24) 

where are the average values of parameters ai, oji, 7 *, respectively, 

and are the bounds of the two virtual state variables, then global 

bounded synehronization is achieved by the network. 

Proof. Let us consider an unweighted complete graph of N HKB oscillators 
interconnected via a full state coupling, that is 


< 


pi 


n—no^—1 ’ 

Pi _ I Tfn ,71 ^ _£L 


> —i") \JU:'i)\ < —T or vice versa 

n—no-—I’l n—no —1 

3 

\J{i,i)\ . 

n—no^ — 1 ’ 

^ n-nn.-l’ I'-'VJ)‘h ^ n-nn.-l 


or Vice versa . 


SI is not true, it is not possible for the virtual 


N 

- c'^aij {xi - Xj), VzG[l,iV] (25) 

i=i 

where Xi G is the state variable of node i and c := since each node 
in a connected complete graph has N — 1 neighbors. The virtual system 
reads 


2 ; = 


2:2 - cA^ 2 :i + cJ 2 f=i Xji 


— azZ 


+ I3z2 - 7 ) 2:2 - a;^ 2 ;i - cNz 2 + c x 




(26) 


where 2 ; G is the state variable of the virtual system and a \= 

^ := T^YA=i^ii 7 := T^Y!i=\lu ^ := T^YA=i^i- The Jacobian of the 
virtual system is: 
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J{t,z) = 


(27) 


-cN 1 

— {2az2Zi + Cj'^) —azl — 3/3z| — cN + 7 

In order to prove global bounded synchronization of the network, we 
need the virtual system to be contracting. In order to do so, we apply the 
algorithm presented in m to Eq. [271 When i = l,j = 2, it is immediate 
to see that statement SI is true, while S2 and S3 are false (c might be in 
general time varying), leading to I J(l, 2)1 > |J(1,1)| and|J(2,l)| < |J(2,2)|. 
When i = 2,j = 1 instead, inequalities to satisfy and CFS depend on the 
sign of a and (5. Supposing without loss of generality that > 0 as 
usually done in literature mm, it is immediate to see that an inequality 
corresponding to the fulfilment of SI needs to be added to the list of CFS 
generated by the algorithm (a worst case scenario is —cN + 7 < 0), and that 
both S2 and S3 are again false, leading to the two same conditions. This 
means that the network achieves global bounded synchronization when the 
following system is satisfied: 


< 1 > cA^ 

^ l2aziZ2 + u)‘^\ < \az‘l + ~ 7 + cN\ 


77 


jY ^ AT 

\ 2 aziZ 2 + < \az\ + 


7 + ciV| 


(28) 


Supposing that the dynamics of the virtual system is bounded, meaning 
that 1^1 (f)| < zi^^^,\z 2 it)\ < Z 2 ^^^ Vt > 0, we can consider the following 
worst case scenario 


which leads to 


— <c< — 

AT ^ c ^ at 

< cA^ - 7 


(29) 


+ 7)^+7 

N 


<c< 


1 

N 


(30) 


and, as a consequence, to Eq. [2ll 

So we can conclude that if the coupling strength c fulfills Eq. [JH the het¬ 
erogeneous network of HKB oscillators overlying a complete graph achieves 
bounded synchronization. □ 


Remark 1. Note that when the number of nodes in the network N is really 
high, then 1. This means that global bounded synchronization is 

achieved when: 


-h 0)2 7 < c < 1 


(31) 
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7.2 Lyapunov theory 

Let T> be the set of diagonal matrices, the set of positive definite diagonal 
matrices and T>~ the set of negative dehnite diagonal matrices. Let us now 
dehne QUAD and QUAD Affine vector fields [3]. 

Definition 3. Given n x n matrices P G , Wi € D, the vector field fi is 
said to be QUAD(P, WQ iff 

{z - w)'^P[fi(t, z) - fi(t, re)] < (z- w)'^Wi{z - w) (32) 

for any z,w G M” and for any t > 0. 

Definition 4. Given n x n matrices P G Wi G P, the vector field fi is 
said to be QUAD(P, tLj) Affine iff fi{t,Xi) = hi{t,Xi) + gi(t,Xi) and 

• hi is QUAD(P,Wi); 

• 3 M < oo : \\gi{t, z )\\2 < M, Vz G M”, Vt > 0 

Let us consider a heterogeneous network of N agents interconnected via 
a linear coupling: 


N 

Xi{t) = fi{t,Xi) - jjr'^aijT{xi- Xj), OO (33) 

* 

where L G Note that this is a generalization of the full state cou¬ 

pling previously introduced, which can be obtained by setting F = ffi. As 
reported in [3] in details, in order to prove global bounded synchronization 
of a network of N nonidentical QUAD Affine systems coupled via a linear 
interaction protocol, we need hi{t,Xi) to be QUAD(P, kU) with Wi G P~ 
for all the nodes in the network at any time instant. However, in a hetero¬ 
geneous network of N HKB oscillators with vector fields described by Eq. 
m regardless of the way we dehne hi and gi it is never possible to satisfy the 
following condition 

{z — w)^P[hi{t, z) — hi{t, re)] < (z — w)'^Wi{z — w) (34) 

with dehnite negative matrices Wi. Indeed, the right-hand term is always 
negative, while the left-hand one can be positive for any value of P > 0. On 
the other hand, in order to avoid conditions on the sign of the matrices Wi, 
it is necessary to write the dynamics of the nodes in the following way 

fi{t,Xi) = hi{t,Xi) + gi{t,Xi) Mi = 1,2,..., N (35) 

with hi(t,z) = hj{t,z) = h{t,z) Mi,j G [1, A"], Mz G M"’, and with all the 
terms gi being bounded, at any time instant. 

In particular, in [3] the authors formalize the following theorem. 
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Theorem 2. Let us consider a heterogeneous network of N agents inter¬ 
connected via a linear coupling as described in Eq. [53 Let us suppose that 
fi{t,Xi) = h{t,Xi) + gi{t,Xi) and that: 

1. the network is made up of N QUAD(P,W) Affine systems, with P G 
V+ and W e V; 

2. r is a positive semidefinite diagonal matrix; 

3. ifW is made up of I ^ [0,n] non-negative elements, which without loss 
of generality can be collected in its Ixl upper-left block, then T is made 
up of I positive elements, where I < I < n, which again without loss of 
generality can be collected in its I x I upper-left block; 

4- 3 0 < M < oo : \\gi{t,Xi )\\2 < M Vi = 1,2, ...,N,\/t > 0. 


Then, we can claim that global bounded synchronization is achieved by 
the network. In particular, if we define matrix Ljg- = {IjVij} as 


Wii ■= < 


1 


LLi-j 

'JTi^ 




if i = j 

if i 7 ^ j and {i,j) are neighbors 
otherwise 


(36) 


we can state that 3 0 < c < oo, e > 0 : limt_).oo vit) < e \/c> c, where 


c = min max 
P,W 


\\2{LM(3)PiTiy ) 


(37) 


with Wi,Pi,Ti representing the I x I upper-left block of matrices W,P,T, 
respectively, and where for a given value of c> c 


y/NM\\P\\2 

min - 

P,W - max {Xm (Wi) - 0 X 2 {Ljg- (g) PiTi ), Am (Wn-i)) 


(38) 


with Wn-i representing the {n — 1) x {n — 1) lower-right block of matrix W 
and with the assumption that 0 X 2 {Lj\f (g) PiTi) > Xm (Wi). 


Proof. See [3]. 


□ 


We can thus derive the following corollary. 


Corollary 1. Let us consider a heterogeneous network of N HKB oscilla¬ 
tors interconnected via a linear coupling. Supposing that the topology of the 
network is simple and undirected, and assuming that 7 i = 7 Vi G [l,N], if 
the coupling strength satisfies the inequality 


c> c = 


mm 


max(W(l,l), 7 P( 2 , 2 )) 


H^(i,i),P(i,i),P(2,2)>o A 2 (Ta^) minj=i ,2 {P{j, j)T{j, j)) 


(39) 
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then global bounded synehronization is achieved by the network. In partieu- 
lar, we can claim that 


lim r]{t) < e 


min 

W{l,l),P{l,l),P(2,2),de>0 


y/NM max (P(l, 1), P(2, 2)) 

4 


(40) 


where 


d, :=c\ 2 {Ljyf) min (P(j, j)r(j, j)) - max (1^(1,1), 7^(2, 2)) (41) 

Proof. First of all we need to write the dynamics of each node in the network 
as fi{t,Xi) = h{t,Xi) + gi{t,Xi). This is possible if we suppose that 7 * = 
7 Vi G [1,-/V] and define: 


h{t,Xi) 


0 

IXi., 


9i{t,Xi) 


Xi2 

-{aixj^ + I3ixj^)xi^ -ujfxi^ 


Then we need to verify whether the nodes in the network are QUAD(P, W) 
Affine systems. In particular, this means that we need h to be QUAD(P, W), 
with P G 2?“*" and VF G P. Therefore, if we define P = diag{P{l, 1), P(2, 2)} 
with P(l, 1), P{2, 2) > 0, VF = diag{W (1,1), VF(2, 2)} and h{t, z) = [0 72 : 2 ]'^ Vz G 
we have to satisfy: 


P(2, 2)7(z2 - W2)^ < VF(1, l)(zi - wif + W{2, 2){z2 - W2f (42) 
Hence, if we choose VF(2, 2) = ^P{2, 2), it possible to reduce Eq. 1421 to 

VF(l,l)(zi-rci)^ > 0 (43) 

which is true for any 1F(1,1) > 0. This means that the first hypoth¬ 
esis of Theorem [2] simply reduces to choosing any P G 'D'^ and IF = 
diag{W{l, 1),^P{2,2)} for any 1F(1,1) > 0. 

Since W G is made up of 2 non-negative elements, we have that 

1 = 1 = 2. Therefore, in order to satisfy the second and the third hypotheses 
of Theorem [21 we need T to be a diagonal positive definite matrix, that is 
r G 'D'^ (note that this is true when the nodes are connected through a full 
state coupling, since it corresponds to T = 12 )- 

Finally, the last hypothesis to satisfy is related to the boundedness of 
the terms gi at any time instant. As already shown before, we have chosen: 


9iit,Xi) 


Xi2 

-{aixj^ + f3ixjjxi2 - iofxi^ 


(44) 


26 









Since the dynamics of each HKB oscillator is bounded m, we can dehne 

Pima. := sup(|xii(t)|), := sup(|xi2(t)|) 

t>0 t>0 

and we define as well pM ■= max* vm ■= max* om ■= 

maxi(|aj|), f3M ■= maxj(|/3j|), com '■= niaxj (Icjjl). Therefore, from Eq. 33] 
we get: 


llffilb < kial + \{<^iXi^ + I3ix‘i^)xi^ 

< Ixjjl + +^i\xii\ 


< (1 + I«*IpL. + \l3i\vLaMma.+^bima. ■= 

Besides, we have that 


(45) 


Mi < (1 + oimPm + ^mx‘m)vm + ^mPm ■— M (46) 

This means that the fourth hypothesis of Theorem [2] is always satished 
in the case of HKB oscillators, and the bound M is dehned in Eq. HH] 

In order to find an easier expression of the minimum value required for 
the coupling strength and of the upper-bound for the error norm, we can 
take advantage of the particular form of matrices P and W: 


P = P2 


P(l,l) 0 

0 E(2,2)J’ 


P(1,1),P(2,2) >0 


W = W2 


'IT(1,1) 0 

0 lP{2,2)\^ 


IT(1,1) > 0 


Therefore, from Eq. |37|we have that the minimum value c for the cou¬ 
pling strength that guarantees bounded synchronization of the network is 
given by Eq. 15^1 while for a given c > c, from Eq. |35| we have that the 
upper bound of the error norm is given by Eq. HOj 

So we can conclude that if c > c > 0, where c is dehned in Eq. [39l global 
bounded synchronization is achieved. □ 


7.3 Numerical validation 

As previously shown for a connected simple undirected heterogeneous net¬ 
work of N HKB oscillators, by making use of contraction theory it is possi¬ 
ble to guarantee bounded synchronization if we suppose that the underlying 
topology is an unweighted complete graph (all-to-all network). On the other 
hand, by making use of Lyapunov theory, bounded synchronization can be 
achieved regardless of the topology and the weights of the interconnections. 
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Figure 9 Underlying topology - simple connected weighted graph 


although an assumption has to be made on one of the parameters of the 
nodes (7* = 7 Vi € [1, A^]). In order to be able to study the most possible 
general case, we consider a weighted random graph of = 5 HKB nonlin¬ 
ear oscillators characterized by 7* = 7 Vi E [1;A^]- In particular, in such 
a random graph the odds of an edge connecting two nodes is worth 60% 
and its weight is randomly picked between 0 and 2 (see Fig. ED- As for the 
parameters and the initial conditions of the nodes, see Table [6j Moreover, 
we set T = 200s. 

Table 6 Numerical simulations - parameters and initial values for a network 
of A^ = 5 HKB oscillators 


Nodes 

ai 

A 

li 

Ui 

Xj(0) 

Node 1 

0.46 

1.16 

0.58 

0.16 

[-1.4,+0.3] 

Node 2 

0.37 

1.20 

0.58 

0.26 

[+1.0,+0.2] 

Node 3 

0.34 

1.73 

0.58 

0.18 

[-1.8,-0.3] 

Node 4 

0.17 

0.31 

0.58 

0.21 

[+0.2,-0.2] 

Node 5 

0.76 

0.76 

0.58 

0.27 

[+1.5,+0.1] 


This scenario leads to pM = 2.6, vm = 0.96, M = 7.6, A 2 {L^f) = 0.4112, 
P(l,l) = 0.077, P(2,2) = 0.077, IU(1,1) = 0.001, IU(2,2) = 0.045 and 
c = 1.4211. Fig. flUl shows xi{t) for all the nodes in the network (blue line: 
node 1, green line: node 2, red line: node 3, cyan line: node 4, magenta line: 
node 5) when connected through a full state coupling protocol with c = 1.45: 
as we can see, our analytical results are confirmed and synchronization is 
achieved by the network. 
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(a) first state variable (b) first state variable - zoom 

Figure 10 First state variable Xjj in a simple connected weighted hetero¬ 
geneous network of N = 5 HKB oscillators - full state coupling (c = 1.45) 


On the other hand, in Fig. [TT] we show that bounded synchronization 
can actually be achieved for smaller values of the coupling strength when 
considering a full state coupling (c = 0.07), and that it can be achieved 
also with the two other coupling protocols presented earlier in this paper 
(ci = C 2 = 0.1 for the partial state coupling and a = b = —1, c = 0.1 for the 
HKB coupling, respectively). Indeed, with this choice of coupling strength, 
the error norm r](t) is roughly bounded by e ~ 2. In Fig. [12] we then show 
the trend of the group synchronization obtained respectively in all the three 
cases: as we can see, after an initial transient Pg{t) reaches a much higher 
value, confirming what observed in [25]. In particular, the trend obtained 
when considering an HKB coupling resembles the most the one obtained in 
the real experiments involving human people: indeed, in this case the group 
synchronization presents persistent oscillations with the highest amplitude, 
as observed in the rocking chairs experiments. 

8 Conclusion 

We have proposed a mathematical model for movement synchronization of 
a group of three or more people. In particular we have considered heteroge¬ 
neous networks of HKB nonlinear oscillators, in which each equation is char¬ 
acterized by a different set of parameters to account for human-to-human 
variability. Three different coupling models, both linear and nonlinear, have 
been investigated, and the effects of adding an external entrainment signal 
have been analyzed. We have found analytical conditions for a connected 
simple undirected network to achieve bounded synchronization when consid¬ 
ering a full state coupling as interaction protocol among the nodes, while we 
have numerically shown that bounded synchronization can be achieved also 
when considering a partial state coupling or a HKB coupling. In particular, 
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Figure 11 Error norm in a simple connected weighted heterogeneous net¬ 
work of = 5 HKB oscillators - magenta solid line; full state coupling 
(c = 0.07), blue solid line: partial state coupling (ci = C 2 = 0.1), black 
dashed line: HKB coupling {a = b = —l,c = 0.1) 




(a) group synchronization 


(b) group synchronization - zoom 


Figure 12 Group synchronization in a simple connected weighted hetero¬ 
geneous network oi N = 5 HKB oscillators - magenta solid line; full state 
coupling (c = 0.07), blue solid line: partial state coupling (ci = C 2 = 0.1), 
black dashed line: HKB coupling {a = b = —1, c = 0.1) 
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we have observed that it is possible to replicate some of the synchroniza¬ 
tion features obtained in the rocking chairs experiments with all the three 
coupling protocols proposed in this paper, although the most realistic one is 
achieved when connecting the nodes through a nonlinear HKB coupling. In¬ 
deed, in this case the group synchronization presents persistent oscillations 
with the highest amplitude, as observed in [25]. 

Some viable extensions of this work include performing real experiments 
involving hand synchronization of more than two players and then choosing 
the coupling protocol that can best capture the observations coming from 
such experiments and explain the onset and features of movement coordi¬ 
nation among the human players. Finally, we wish to prove global bounded 
synchronization with any kind of coupling protocol, both linear and nonlin¬ 
ear, also in the more general case of directed networks. 
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